% $Id: kalman.tex 204 2009-04-20 01:59:36Z jsibert $
%
% Author: David Fournier
% Copyright (c) 2008 Regents of the University of California
%
%\input docmacpdf.tex
%\def\aone{\hbox{$a$\kern -1pt $1$}}
%\def\Pone{\hbox{$P$\kern -1pt $1$}}
%\chapno=1
\def\aone{a_{t|t-1}}
\def\Pone{P_{t|t-1}}
\mysection{The Kalman Filter}

The Kalman filter is a device for estimating parameters in
a class of ``time-series'' like models which are put into state-space
form. We have used the 
notation from Harvey, chapter 3.
\X{state space form}
\X{Kalman filter}
The general state space form is an  multivariate time series
$$y_t=Z_t\alpha_t+d_t+\epsilon_t$$ 
where $Z_t$ is an $N\times m$ matrix, $d_t$ is an $N$ dimensional
vector, $y_t$ is an $N$ dimensional vector and $\epsilon_t$ is
a set of serially uncorrelated $N$ dimensional random vectors 
with mean $0$ and correlation $H_t$.
The elements of $\alpha_t$ are not observable but are
assumed to be generated by a first order Markov process
$$\alpha_t=T_t\alpha_{t-1}+c_t+R_t\eta_t$$
where $T_t$ is an $m\times m$ matrix, $c_t$ is an $m\times 1$
vector, $R_t$ is an $m\times g$ matrix and $\eta_t$
is a $g\times 1$ vector of serially uncorrelated 
random vectors with mean $0$ and covariance matrix $H_t$.
The specification of the state space system is completed by
two further assumptions:
\smallskip
1. The initial state vector $\alpha_0$ has a mean of $a_0$ and a 
variance of $P_0$.
\smallskip
2. The random vectors $\epsilon_t$ and $\eta_t$ are uncorrelated
with each other and uncorrelated with the initial state.

In applications of the model many of the
parameters $Z_t$, $d_t$, $H_t$, $T_t$, $c_t$, $R_t$, and $Q_t$
may be independent of $t$ in which case we will write them
without the subscript. Also $R$ may be the identity matrix
in which case we will omit it.

As a simple example of such a model consider the (two dimensional)
random walk observed with error example considered below.
%$$ \eqalign{
 \begin{eqnarray}
   \alpha_t&=&\alpha_{t-1}+\eta_t\cr
   y_t&=&\alpha_t+\epsilon_t\cr
 \end{eqnarray}
%} $$
For this model  the following parameters are fixed.
$$ T=\left(\begin{matrix}1 & 0\\
                    0 & 1 \end{matrix}\right) 
  \quad
   Z=\left(\begin{matrix}1 & 0\cr
                    0 & 1 \end{matrix}\right) 
  \quad
   R=\left(\begin{matrix}1 & 0\cr
                    0 & 1\end{matrix}\right) 
  \quad
  d=(0,0)
  \quad
  c=(0,0)
$$
while the covariance matrices $Q$ and $H$ are estimated. Their
true values used for simulating the data were. The initial
value of $a$ is $(0,0)$.
$$
   Q=\left(\begin{matrix}1 & 0.8\cr
                   0.8 & 1 \end{matrix}\right)
  \quad
   H=\left(\begin{matrix}3 & -2.5\cr
                   -2.5 & 3 \end{matrix}\right)
$$
 
\mysection{Equations for the Kalman filter}
For a moment go back to the general state space model.
Given $a_0$, $P_0$ we recursively calculate the a number of quantities
via the relationships
%$$\eqalign{
 \begin{eqnarray}
      \aone&=&T_ta_{t-1}+c_t\cr
      \Pone&=&T_tP_{t-1}T^\prime_t +R_tQ_tR^\prime_t\cr
      v_t&=&y_t-Z_t\aone-d_t\cr
      F_t&=&Z_t\Pone Z_t^\prime+H_t\cr
      a_t&=&\aone-\Pone Z_t^\prime F_t^{-1}v_t\cr
      P_t&=&\Pone-\Pone Z_t^\prime F_t^{-1}Z_t\Pone
 \end{eqnarray}
%} $$
The log likelihood function for the models parameters is
given by
$$\log L = -{NT\over 2}\log 2\pi -0.5 \sum_{t=1}^T \log|F_t|
        -0.5 \sum_{t=1}^T v_t F_t^{-1} v_t$$

The {\tt TPL} file for the random walk kalman filter code
follows
\beginexample
DATA_SECTION
  init_int nt
  init_int N
  init_int m
  int m1
  init_matrix Y(1,nt,1,N)
  matrix P0(1,m,1,m)
 !! P0.initialize(); 
 !! m1=m*(m+1)/2; 
PARAMETER_SECTION
  init_bounded_vector Qcoff(1,m1,-10.,10.1)
  init_bounded_vector Hcoff(1,m1,-10.,10.1)
  init_vector a0(1,m)
  matrix T(1,m,1,m)
  matrix TT(1,m,1,m)
  vector d(1,N)
  vector c(1,m)
  matrix chQ(1,m,1,m)
  sdreport_matrix Q(1,m,1,m)
  matrix chH(1,N,1,N)
  sdreport_matrix H(1,N,1,N)
  matrix Z(1,N,1,m)
  matrix TZ(1,m,1,N)
  objective_function_value f
 LOCAL_CALCS
   d.initialize();
   c.initialize();
   Z.initialize();
   Z(1,1)=1; Z(2,2)=1;
   T.initialize();
   T(1,1)=1; T(2,2)=1;
   TZ=trans(Z);
   TT=trans(T);
PROCEDURE_SECTION
  setup_Q();
  setup_H();
  f+=kalman_filter();
  cout << " f = " << f << endl;

FUNCTION setup_Q
  chQ.initialize();
  int ii=1;
  for (int i=1;i<=m;i++)
    for (int j=1;j<=i;j++)
      chQ(i,j)=Qcoff(ii++);
  Q=chQ*trans(chQ);    
FUNCTION setup_H
  chH.initialize();
  int ii=1;
  for (int i=1;i<=N;i++)
    for (int j=1;j<=i;j++)
      chH(i,j)=Hcoff(ii++);
  H=chH*trans(chH);    

FUNCTION dvariable kalman_filter(void)  
  dvar3_array P(0,nt,1,m,1,m);
  dvar3_array P1(1,nt,1,m,1,m);
  dvar3_array F(1,nt,1,N,1,N);
  dvar3_array Finv(1,nt,1,N,1,N);
  dvar_matrix Ptemp(1,m,1,m);
  dvar_matrix a(0,nt,1,m);
  dvar_matrix a1(1,nt,1,m);
  dvar_matrix v(1,nt,1,N);
  a(0)=a0;
  P(0)=P0;
  // This is the Kalman filter recursion. The objects tmp1
  // and tmp2 hold common calculations to optimize a bit
  int t;
  for (t=1;t<=nt;t++)
  {
    a1(t)=T*a(t-1)+c;
    P1(t)=T*P(t-1)*TT+Q;
    dvar_vector pred_y=Z*a1(t)+d;
    v(t)=Y(t)-pred_y;
    dvar_matrix tmp1=P1(t)*TZ;
    F(t)=Z*tmp1+H;
    Finv(t)=inv(F(t));
    dvar_matrix tmp2= tmp1*Finv(t);
    P(t)=P1(t)-tmp2*Z*P1(t);
    a(t)=a1(t)+tmp2*v(t);
  }
  int sgn=0;
  dvariable f=0.0;
  for (t=1;t<=nt;t++)
    f+=ln_det(F(t),sgn)+v(t)*Finv(t)*v(t);
  return f;  
TOP_OF_MAIN_SECTION
  arrmblsize=20000000;
  gradient_structure::set_CMPDIF_BUFFER_SIZE(3000000);
  gradient_structure::set_GRADSTACK_BUFFER_SIZE(1000000);
 
\endexample 
This example was deliberatley not optimized as much as it could be
in order to retain the flavour of the more general state space
problem. For example since $T$ is the identity matrix and $c$ is the
zero vector the line of code 
\beginexample
    a1(t)=T*a(t-1)+c;
\endexample 
reduces to
\beginexample
    a1(t)=a(t-1);
\endexample 
The parameters being estimated are  $a_0$, $Q$, and $H$. 

To parameterize the covariance matrices the Choleski decomposition
parameterization was used. This ensures that the covariance matrices are
positive (semi) definite. The technique can be seen in the function
{\tt setup\_Q}. 
The lower triangular matrix {\tt ch\_Q} is filled with parameters
from a bounded vector
\beginexample
FUNCTION setup_Q
  chQ.initialize();
  int ii=1;
  for (int i=1;i<=m;i++)
    for (int j=1;j<=i;j++)
      chQ(i,j)=Qcoff(ii++);
  Q=chQ*trans(chQ);   // chQ is the choleski decomposition of Q 
\endexample
Notice that the bounded vector {\tt Qcoff} has slightly asymmetric
bounds. This is a simple way to ensure that its initial value is not
identically 0 which would lead to a singular covariance matrix.
\beginexample
  init_bounded_vector Qcoff(1,m1,-10.,10.1)
\endexample

The model parameters, standard deviations and 
corellations are reproduced from the
standadard ADMB report.

\beginexample
\obeylines
\stt
 index   name    value      std dev    7      8      9     10     11     12     13     14     15      16   
    7   a0    -1.1682e+00 9.0191e-01  1.000
    8   a0     1.2218e+00 8.6442e-01  0.352  1.000
    9   Q      9.9468e-01 1.0862e-01  0.059 -0.006  1.000
   10   Q      7.8808e-01 7.8737e-02  0.038  0.028  0.683  1.000
   11   Q      7.8808e-01 7.8737e-02  0.038  0.028  0.683  1.000  1.000
   12   Q      8.7279e-01 9.6118e-02 -0.018  0.069  0.185  0.721  0.721  1.000
   13   H      3.1352e+00 1.8123e-01 -0.015 -0.007 -0.305 -0.136 -0.136 -0.018  1.000
   14   H     -2.7119e+00 1.4922e-01 -0.021  0.001 -0.102 -0.238 -0.238 -0.139 -0.692  1.000
   15   H     -2.7119e+00 1.4922e-01 -0.021  0.001 -0.102 -0.238 -0.238 -0.139 -0.692  1.000  1.000
   16   H      3.2264e+00 1.7936e-01  0.015 -0.029 -0.031 -0.121 -0.121 -0.249  0.370 -0.698 -0.698  1.000
\endexample
\mysection{Parameterizing the covariance matrix parameterizations}

The Choleski decomposition parameterization merely ensures that the 
matrix is positive semi-definite. By adding a small positive number to the
diagonal elements one can ensure that the covariance matrix is
positive definite and can speed up and improve the stability of the
estimation. Of course what is meant by small will depend on the
particular problem being considered. 
A modifed form of the routine {\tt setup\_Q}
follows.
\beginexample
FUNCTION setup_Q
  int i;
  chQ.initialize();
  int ii=1;
  for (i=1;i<=m;i++)
    for (int j=1;j<=i;j++)
      chQ(i,j)=Qcoff(ii++);
  Q=chQ*trans(chQ);   // chQ is the choleski decomposition of Q 
  for (i=1;i<=m;i++)
    Q(i,i)+=0.1;   // make Q positive definite
\endexample
Performing this modification for the present model for both $Q$ and
$H$ causes the program to converge about twice as fast.

